knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE 

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_SensitivityBarriersFacilitators_', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan', # todo: do we need this??
    'barriers_self_cw', 
    'facilitators_self_cw', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'barriers_self_cb', 
    'facilitators_self_cb', 
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily individual's experienced persuasion",  
    "Daily partner's experienced persuasion", 
    "Daily individual's experienced pressure", 
    "Daily partner's experienced pressure", 
    "Daily individual's experienced pushing", 
    "Daily partner's experienced pushing", 
    "Day", 
    "Own action plan",
    'Partner action plan',
    "Daily barriers",
    "Daily facilitators",
    "Daily wear time",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean individual's experienced persuasion", 
    "Mean partner's experienced persuasion", 
    "Mean individual's experienced pressure", 
    "Mean partner's experienced pressure", 
    "Mean individual's experienced pushing", 
    "Mean partner's experienced pushing", 
    'Mean barriers',
    "Mean facilitators",
    "Mean wear time"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily individual's experienced persuasion)", 
  "sd(Daily partner's experienced persuasion)", # OR partner received
  "sd(Daily individual's experienced pressure)", 
  "sd(Daily partner's experienced pressure)", 
  "sd(Daily individual's experienced pushing)", 
  "sd(Daily partner's experienced pushing)", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,13),
  "Between-Person Effects" = c(14,22),
  "Random Effects" = c(23, 29), 
  "Additional Parameters" = c(30,30)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,13+5),
  "Between-Person Effects" = c(14+5,22+5),
  "Random Effects" = c(23+5, 29+5), 
  "Additional Parameters" = c(30+5,30+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
## Start sampling
pa_sub_digest <- digest::digest(pa_sub)
if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.35 [1.30, 1.40]         1.16      0.74
##  persuasion_partner_cw 1.08 [1.05, 1.12]         1.04      0.93
##       pressure_self_cw 1.05 [1.03, 1.09]         1.02      0.95
##    pressure_partner_cw 1.07 [1.04, 1.11]         1.03      0.94
##        pushing_self_cw 1.03 [1.01, 1.08]         1.02      0.97
##     pushing_partner_cw 1.14 [1.11, 1.19]         1.07      0.87
##     persuasion_self_cb 1.05 [1.02, 1.09]         1.02      0.95
##  persuasion_partner_cb 3.74 [3.56, 3.93]         1.93      0.27
##       pressure_self_cb 3.66 [3.48, 3.85]         1.91      0.27
##    pressure_partner_cb 2.37 [2.27, 2.48]         1.54      0.42
##        pushing_self_cb 2.32 [2.22, 2.43]         1.52      0.43
##     pushing_partner_cb 3.24 [3.09, 3.40]         1.80      0.31
##       barriers_self_cw 3.26 [3.10, 3.42]         1.80      0.31
##       barriers_self_cb 1.25 [1.21, 1.30]         1.12      0.80
##   facilitators_self_cw 1.10 [1.07, 1.14]         1.05      0.91
##   facilitators_self_cb 1.26 [1.22, 1.31]         1.12      0.79
##                    day 1.08 [1.05, 1.12]         1.04      0.92
##  Tolerance 95% CI
##      [0.71, 0.77]
##      [0.89, 0.95]
##      [0.91, 0.97]
##      [0.90, 0.96]
##      [0.92, 0.99]
##      [0.84, 0.90]
##      [0.92, 0.98]
##      [0.25, 0.28]
##      [0.26, 0.29]
##      [0.40, 0.44]
##      [0.41, 0.45]
##      [0.29, 0.32]
##      [0.29, 0.32]
##      [0.77, 0.82]
##      [0.87, 0.93]
##      [0.76, 0.82]
##      [0.89, 0.95]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         91%
##     lognormal          9%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         84%
##               beta-binomial          9%
##                   lognormal          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 7 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 4, observations = 1776, p-value =
## 0.72
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.003378378
## sample estimates:
## outlier frequency (expected: 0.00166103603603604 ) 
##                                        0.002252252
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report,
## rope_range = rope_range_continuous, : Coefficients were exponentiated.
## Double check if this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 2.50*** 0.62 [ 1.51, 4.14] 1.000 [0.84, 1.20] 0.003 26.868 Strong Evidence 1.000 30408 29858 49.95*** 3.90 [42.68, 58.63] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 11207 20136
Within-Person Effects
Daily individual’s experienced persuasion 1.46*** 0.14 [ 1.23, 1.79] 1.000 [0.84, 1.20] 0.012 >100 Overwhelming Evidence 1.000 48576 28761 1.02 0.03 [ 0.97, 1.09] 0.811 [0.93, 1.08] 0.960 0.016 Very Strong Evidence for Null 1.000 24252 30545
Daily partner’s experienced persuasion 1.32** 0.13 [ 1.10, 1.64] 0.999 [0.84, 1.20] 0.141 4.547 Moderate Evidence 1.000 42058 28118 1.02 0.02 [ 0.98, 1.07] 0.850 [0.93, 1.08] 0.984 0.015 Very Strong Evidence for Null 1.000 38268 30627
Daily individual’s experienced pressure 0.80 0.20 [ 0.47, 1.49] 0.805 [0.84, 1.20] 0.353 0.206 Moderate Evidence for Null 1.000 31438 21163 0.88 0.06 [ 0.77, 1.01] 0.970 [0.93, 1.08] 0.216 0.073 Strong Evidence for Null 1.000 34797 28176
Daily partner’s experienced pressure 1.52 0.60 [ 0.82, 5.09] 0.898 [0.84, 1.20] 0.222 0.314 Weak Evidence for Null 1.000 24946 22538 0.95 0.04 [ 0.86, 1.04] 0.864 [0.93, 1.08] 0.693 0.015 Very Strong Evidence for Null 1.000 44931 32944
Daily individual’s experienced pushing 1.17 0.17 [ 0.88, 1.60] 0.867 [0.84, 1.20] 0.553 0.131 Moderate Evidence for Null 1.000 36871 28172 0.97 0.03 [ 0.91, 1.03] 0.808 [0.93, 1.08] 0.940 0.012 Very Strong Evidence for Null 1.000 38610 36953
Daily partner’s experienced pushing 1.60** 0.27 [ 1.18, 2.47] 0.998 [0.84, 1.20] 0.031 6.066 Moderate Evidence 1.000 28004 22866 0.97 0.03 [ 0.91, 1.02] 0.875 [0.93, 1.08] 0.920 0.016 Very Strong Evidence for Null 1.000 46513 33550
Day 0.97 0.22 [ 0.63, 1.50] 0.550 [0.84, 1.20] 0.579 0.110 Moderate Evidence for Null 1.000 72916 28349 1.00 0.07 [ 0.88, 1.14] 0.501 [0.93, 1.08] 0.737 0.008 Very Strong Evidence for Null 1.000 69173 29875
Own action plan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Partner action plan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.71*** 0.07 [ 0.58, 0.86] 1.000 [0.84, 1.20] 0.054 12.137 Strong Evidence 1.000 63912 33768 1.22*** 0.05 [ 1.13, 1.31] 1.000 [0.93, 1.08] 0.001 >100 Overwhelming Evidence 1.000 66926 31862
Daily facilitators 1.31** 0.14 [ 1.07, 1.61] 0.995 [0.84, 1.20] 0.190 1.631 Weak Evidence 1.000 66058 31798 1.23*** 0.03 [ 1.17, 1.29] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 77685 29208
Daily wear time NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.83 0.45 [ 0.28, 2.44] 0.636 [0.84, 1.20] 0.242 0.291 Moderate Evidence for Null 1.000 19867 26740 1.06 0.17 [ 0.77, 1.45] 0.647 [0.93, 1.08] 0.345 0.016 Very Strong Evidence for Null 1.000 20717 27262
Mean partner’s experienced persuasion 1.27 0.69 [ 0.43, 3.69] 0.669 [0.84, 1.20] 0.235 0.296 Moderate Evidence for Null 1.000 20034 26527 0.92 0.14 [ 0.67, 1.26] 0.698 [0.93, 1.08] 0.330 0.017 Very Strong Evidence for Null 1.000 19148 27680
Mean individual’s experienced pressure 0.32 0.24 [ 0.07, 1.28] 0.945 [0.84, 1.20] 0.059 1.239 Weak Evidence 1.000 29411 31320 0.65 0.25 [ 0.30, 1.42] 0.865 [0.93, 1.08] 0.082 0.044 Strong Evidence for Null 1.000 8942 16260
Mean partner’s experienced pressure 0.25* 0.18 [ 0.05, 1.00] 0.975 [0.84, 1.20] 0.032 2.377 Weak Evidence 1.000 29140 30083 0.49 0.19 [ 0.23, 1.09] 0.960 [0.93, 1.08] 0.032 0.121 Moderate Evidence for Null 1.000 8685 15802
Mean individual’s experienced pushing 1.09 0.87 [ 0.22, 5.34] 0.544 [0.84, 1.20] 0.177 0.406 Weak Evidence for Null 1.000 25331 28224 1.48 0.42 [ 0.83, 2.64] 0.913 [0.93, 1.08] 0.082 0.045 Strong Evidence for Null 1.000 11775 20612
Mean partner’s experienced pushing 1.48 1.20 [ 0.30, 7.25] 0.684 [0.84, 1.20] 0.158 0.456 Weak Evidence for Null 1.000 25349 29917 1.70 0.49 [ 0.95, 3.05] 0.963 [0.93, 1.08] 0.041 0.086 Strong Evidence for Null 1.000 11460 19048
Mean barriers 0.54** 0.13 [ 0.34, 0.86] 0.996 [0.84, 1.20] 0.032 3.859 Moderate Evidence 1.000 75540 30068 0.94 0.07 [ 0.80, 1.09] 0.802 [0.93, 1.08] 0.517 0.014 Very Strong Evidence for Null 1.000 55221 34199
Mean facilitators 1.28 0.22 [ 0.91, 1.79] 0.921 [0.84, 1.20] 0.349 0.231 Moderate Evidence for Null 1.000 59512 31381 1.10 0.06 [ 0.99, 1.21] 0.963 [0.93, 1.08] 0.368 0.061 Strong Evidence for Null 1.000 35442 32659
Mean wear time NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.23 0.19 [0.92, 1.67] NA NA NA NA NA 1.000 22060 28077 0.29 0.04 [0.22, 0.40] NA NA NA NA NA 1.000 16410 23783
sd(Daily individual’s experienced persuasion) 0.16 0.12 [0.01, 0.43] NA NA NA NA NA 1.000 17036 20110 0.11 0.03 [0.07, 0.17] NA NA NA NA NA 1.000 22724 25370
sd(Daily partner’s experienced persuasion) 0.22 0.14 [0.01, 0.52] NA NA NA NA NA 1.000 14165 17113 0.06 0.03 [0.01, 0.12] NA NA NA NA NA 1.000 15863 12381
sd(Daily individual’s experienced pressure) 0.44 0.41 [0.02, 1.84] NA NA NA NA NA 1.000 12861 21861 0.09 0.08 [0.00, 0.33] NA NA NA NA NA 1.000 15519 20700
sd(Daily partner’s experienced pressure) 0.82 0.62 [0.05, 2.38] NA NA NA NA NA 1.000 12925 16955 0.05 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 22490 20987
sd(Daily individual’s experienced pushing) 0.42 0.23 [0.04, 0.94] NA NA NA NA NA 1.000 10965 12664 0.06 0.04 [0.00, 0.15] NA NA NA NA NA 1.000 13431 16351
sd(Daily partner’s experienced pushing) 0.42 0.29 [0.03, 1.15] NA NA NA NA NA 1.000 10306 17735 0.05 0.04 [0.00, 0.13] NA NA NA NA NA 1.000 13963 16152
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.66 0.01 [0.63, 0.68] NA NA NA NA NA 1.000 64050 28796 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.76), b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.71), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.8). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Hurdle part of the model on the left, non-zero part towards the right side of the table

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Scanning ttf files in C:\Windows/Fonts, C:\Users\pascku\AppData\Local/Microsoft/Windows/Fonts ...
## Extracting .afm files from .ttf files...
## C:\Windows\Fonts\Candara.ttf : Candara already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarab.ttf : Candara-Bold already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarai.ttf : Candara-Italic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaral.ttf : Candara-Light already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarali.ttf : Candara-LightItalic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaraz.ttf : Candara-BoldItalic already registered in fonts database. Skipping.
## Found FontName for 0 fonts.
## Scanning afm files in C:/Users/pascku/AppData/Local/R/cache/R/renv/cache/v5/windows/R-4.4/x86_64-w64-mingw32/extrafontdb/1.0/a861555ddec7451c653b40e713166c6f/extrafontdb/metrics
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 134 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 208 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.00764

$persuasion_partner_cw

## Warning: Removed 108 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 203 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.00698

$pressure_self_cw

## Warning: Removed 62 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 90 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.0126

$pressure_partner_cw

## Warning: Removed 8 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.029

$pushing_self_cw

## Warning: Removed 4 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.00901

$pushing_partner_cw

## Picking joint bandwidth of 0.0129

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Comparing effect size of pressure and pushing

hypothesis(pa_sub, "pressure_self_cw < pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper
## 1 (pressure_self_cw... < 0     -0.1      0.08    -0.23     0.03
##   Evid.Ratio Post.Prob Star
## 1       8.89       0.9     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner + 
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.12 [1.08, 1.18]         1.06      0.89
##  persuasion_partner_cw 1.16 [1.11, 1.22]         1.08      0.86
##       pressure_self_cw 1.07 [1.04, 1.14]         1.04      0.93
##    pressure_partner_cw 1.10 [1.07, 1.16]         1.05      0.91
##        pushing_self_cw 1.17 [1.13, 1.23]         1.08      0.85
##     pushing_partner_cw 1.17 [1.12, 1.23]         1.08      0.86
##       pressure_self_cb 3.47 [3.25, 3.71]         1.86      0.29
##    pressure_partner_cb 3.62 [3.39, 3.88]         1.90      0.28
##       barriers_self_cw 1.31 [1.25, 1.38]         1.14      0.76
##       barriers_self_cb 1.13 [1.09, 1.19]         1.06      0.89
##   facilitators_self_cw 1.30 [1.25, 1.37]         1.14      0.77
##   facilitators_self_cb 1.11 [1.07, 1.17]         1.06      0.90
##                    day 1.06 [1.03, 1.13]         1.03      0.94
##       weartime_self_cw 1.06 [1.03, 1.12]         1.03      0.95
##       weartime_self_cb 1.21 [1.16, 1.28]         1.10      0.82
##  Tolerance 95% CI
##      [0.85, 0.92]
##      [0.82, 0.90]
##      [0.88, 0.96]
##      [0.86, 0.94]
##      [0.81, 0.89]
##      [0.82, 0.89]
##      [0.27, 0.31]
##      [0.26, 0.30]
##      [0.72, 0.80]
##      [0.84, 0.92]
##      [0.73, 0.80]
##      [0.85, 0.93]
##      [0.89, 0.97]
##      [0.89, 0.97]
##      [0.78, 0.86]
## 
## Moderate Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 6.52 [6.07, 7.01]         2.55      0.15
##  persuasion_partner_cb 7.00 [6.51, 7.53]         2.65      0.14
##        pushing_self_cb 5.63 [5.25, 6.05]         2.37      0.18
##     pushing_partner_cb 5.53 [5.15, 5.95]         2.35      0.18
##  Tolerance 95% CI
##      [0.14, 0.16]
##      [0.13, 0.15]
##      [0.17, 0.19]
##      [0.17, 0.19]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 15, observations = 1594, p-value =
## 0.02
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0006273526 0.0062735257
## sample estimates:
## outlier frequency (expected: 0.003099121706399 ) 
##                                      0.009410289
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report,
## : Coefficients were exponentiated. Double check if this was intended.
summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 122.68*** 7.54 [108.65, 138.51] 1.000 [0.94, 1.06] 0.000 >100 Overwhelming Evidence 1.000 9035 17011
Within-Person Effects
Daily individual’s experienced persuasion 1.02 0.02 [ 0.98, 1.06] 0.865 [0.94, 1.06] 0.990 0.013 Very Strong Evidence for Null 1.000 31618 30796
Daily partner’s experienced persuasion 1.03 0.02 [ 0.99, 1.07] 0.946 [0.94, 1.06] 0.966 0.026 Very Strong Evidence for Null 1.000 40169 30447
Daily individual’s experienced pressure 0.95 0.04 [ 0.86, 1.05] 0.857 [0.94, 1.06] 0.583 0.013 Very Strong Evidence for Null 1.000 41124 28141
Daily partner’s experienced pressure 0.96 0.04 [ 0.89, 1.04] 0.833 [0.94, 1.06] 0.744 0.010 Very Strong Evidence for Null 1.000 52856 30937
Daily individual’s experienced pushing 1.03 0.03 [ 0.97, 1.09] 0.820 [0.94, 1.06] 0.901 0.012 Very Strong Evidence for Null 1.000 40134 32844
Daily partner’s experienced pushing 1.01 0.02 [ 0.96, 1.06] 0.635 [0.94, 1.06] 0.985 0.007 Very Strong Evidence for Null 1.000 51637 30560
Day 0.97 0.05 [ 0.88, 1.07] 0.746 [0.94, 1.06] 0.706 0.007 Very Strong Evidence for Null 1.000 75462 28919
Own action plan NA NA NA NA NA NA NA NA NA NA NA
Partner action plan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.99 0.02 [ 0.94, 1.03] 0.717 [0.94, 1.06] 0.977 0.006 Very Strong Evidence for Null 1.000 74133 29983
Daily facilitators 1.06** 0.02 [ 1.02, 1.10] 0.999 [0.94, 1.06] 0.532 0.570 Weak Evidence for Null 1.000 75238 29620
Daily wear time 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.06] 1.000 1.397 Weak Evidence 1.000 77933 30606
Between-Person Effects
Mean individual’s experienced persuasion 1.01 0.16 [ 0.73, 1.39] 0.522 [0.94, 1.06] 0.302 0.015 Very Strong Evidence for Null 1.001 7858 14370
Mean partner’s experienced persuasion 0.88 0.14 [ 0.63, 1.22] 0.792 [0.94, 1.06] 0.215 0.021 Very Strong Evidence for Null 1.001 7696 13997
Mean individual’s experienced pressure 1.08 0.20 [ 0.75, 1.55] 0.658 [0.94, 1.06] 0.243 0.010 Very Strong Evidence for Null 1.000 11051 18838
Mean partner’s experienced pressure 1.03 0.18 [ 0.72, 1.46] 0.560 [0.94, 1.06] 0.271 0.011 Very Strong Evidence for Null 1.001 9952 17131
Mean individual’s experienced pushing 1.10 0.26 [ 0.69, 1.77] 0.656 [0.94, 1.06] 0.194 0.014 Very Strong Evidence for Null 1.001 10677 19181
Mean partner’s experienced pushing 1.29 0.30 [ 0.81, 2.07] 0.864 [0.94, 1.06] 0.114 0.024 Very Strong Evidence for Null 1.001 10415 17961
Mean barriers 0.79*** 0.04 [ 0.71, 0.88] 1.000 [0.94, 1.06] 0.001 42.237 Very Strong Evidence 1.000 60663 32066
Mean facilitators 0.95 0.04 [ 0.88, 1.03] 0.903 [0.94, 1.06] 0.613 0.022 Very Strong Evidence for Null 1.000 38555 31770
Mean wear time 1.00* 0.00 [ 1.00, 1.00] 0.984 [0.94, 1.06] 1.000 0.094 Strong Evidence for Null 1.000 32032 31012
Random Effects
sd(Intercept) 0.33 0.04 [0.26, 0.44] NA NA NA NA NA 1.000 11399 17633
sd(Daily individual’s experienced persuasion) 0.05 0.02 [0.02, 0.09] NA NA NA NA NA 1.000 16832 13116
sd(Daily partner’s experienced persuasion) 0.05 0.02 [0.00, 0.10] NA NA NA NA NA 1.000 10335 10827
sd(Daily individual’s experienced pressure) 0.06 0.06 [0.00, 0.23] NA NA NA NA NA 1.000 14491 19697
sd(Daily partner’s experienced pressure) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 23179 19565
sd(Daily individual’s experienced pushing) 0.09 0.03 [0.02, 0.16] NA NA NA NA NA 1.000 11366 8683
sd(Daily partner’s experienced pushing) 0.05 0.03 [0.00, 0.12] NA NA NA NA NA 1.000 12918 17397
Additional Parameters
sigma 0.54 0.01 [0.52, 0.56] NA NA NA NA NA 1.000 56096 29752 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.84), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.77). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.18 [1.14,  1.24]         1.09      0.85
##  persuasion_partner_cw 1.14 [1.10,  1.20]         1.07      0.88
##       pressure_self_cw 1.13 [1.09,  1.19]         1.06      0.89
##    pressure_partner_cw 1.08 [1.04,  1.14]         1.04      0.93
##        pushing_self_cw 1.25 [1.20,  1.32]         1.12      0.80
##     pushing_partner_cw 1.15 [1.11,  1.21]         1.07      0.87
##       pressure_self_cb 4.94 [4.62,  5.30]         2.22      0.20
##    pressure_partner_cb 4.94 [4.62,  5.30]         2.22      0.20
##       barriers_self_cw 1.28 [1.23,  1.35]         1.13      0.78
##       barriers_self_cb 1.10 [1.06,  1.16]         1.05      0.91
##   facilitators_self_cw 1.28 [1.22,  1.34]         1.13      0.78
##   facilitators_self_cb 1.08 [1.04,  1.14]         1.04      0.93
##                    day 1.07 [1.04,  1.13]         1.04      0.93
##  Tolerance 95% CI
##      [0.80, 0.88]
##      [0.83, 0.91]
##      [0.84, 0.92]
##      [0.88, 0.96]
##      [0.76, 0.83]
##      [0.83, 0.90]
##      [0.19, 0.22]
##      [0.19, 0.22]
##      [0.74, 0.81]
##      [0.86, 0.94]
##      [0.74, 0.82]
##      [0.88, 0.96]
##      [0.88, 0.96]
## 
## Moderate Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##  persuasion_partner_cb 9.86 [9.17, 10.61]         3.14      0.10
##        pushing_self_cb 7.30 [6.80,  7.84]         2.70      0.14
##     pushing_partner_cb 7.31 [6.81,  7.85]         2.70      0.14
##  Tolerance 95% CI
##      [0.09, 0.11]
##      [0.13, 0.15]
##      [0.13, 0.15]
## 
## High Correlation
## 
##                Term   VIF    VIF 95% CI Increased SE Tolerance
##  persuasion_self_cb 10.02 [9.32, 10.78]         3.17      0.10
##  Tolerance 95% CI
##      [0.09, 0.11]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         47%
##        cauchy         38%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 2 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 10, observations = 1776, p-value =
## 0.003607
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.00270331 0.01033048
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.005630631
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Warning in compute_rope(range = rope_range): Collinearity detected. Some
## VIFs are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF    VIF 95% CI Increased SE Tolerance
persuasion_self_cw 1.18 [1.14,  1.24]         1.09      0.85

persuasion_partner_cw 1.14 [1.10, 1.20] 1.07 0.88 pressure_self_cw 1.13 [1.09, 1.19] 1.06 0.89 pressure_partner_cw 1.08 [1.04, 1.14] 1.04 0.93 pushing_self_cw 1.25 [1.20, 1.32] 1.12 0.80 pushing_partner_cw 1.15 [1.11, 1.21] 1.07 0.87 pressure_self_cb 4.94 [4.62, 5.30] 2.22 0.20 pressure_partner_cb 4.94 [4.62, 5.30] 2.22 0.20 barriers_self_cw 1.28 [1.23, 1.35] 1.13 0.78 barriers_self_cb 1.10 [1.06, 1.16] 1.05 0.91 facilitators_self_cw 1.28 [1.22, 1.34] 1.13 0.78 facilitators_self_cb 1.08 [1.04, 1.14] 1.04 0.93 day 1.07 [1.04, 1.13] 1.04 0.93 Tolerance 95% CI [0.80, 0.88] [0.83, 0.91] [0.84, 0.92] [0.88, 0.96] [0.76, 0.83] [0.83, 0.90] [0.19, 0.22] [0.19, 0.22] [0.74, 0.81] [0.86, 0.94] [0.74, 0.82] [0.88, 0.96] [0.88, 0.96]

Moderate Correlation

              Term  VIF    VIF 95% CI Increased SE Tolerance

persuasion_partner_cb 9.86 [9.17, 10.61] 3.14 0.10 pushing_self_cb 7.30 [6.80, 7.84] 2.70 0.14 pushing_partner_cb 7.31 [6.81, 7.85] 2.70 0.14 Tolerance 95% CI [0.09, 0.11] [0.13, 0.15] [0.13, 0.15]

High Correlation

           Term   VIF    VIF 95% CI Increased SE Tolerance

persuasion_self_cb 10.02 [9.32, 10.78] 3.17 0.10 Tolerance 95% CI [0.09, 0.11]

## Sampling priors, please wait...
summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.80*** 0.11 [ 3.59, 4.02] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 7659 15385
Within-Person Effects
Daily individual’s experienced persuasion 0.01 0.02 [-0.04, 0.06] 0.647 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 56219 31125
Daily partner’s experienced persuasion 0.01 0.02 [-0.04, 0.06] 0.682 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 59605 30043
Daily individual’s experienced pressure -0.05 0.06 [-0.18, 0.08] 0.768 [-0.11, 0.11] 0.831 0.006 Very Strong Evidence for Null 1.000 48822 27686
Daily partner’s experienced pressure -0.06 0.06 [-0.19, 0.06] 0.839 [-0.11, 0.11] 0.797 0.009 Very Strong Evidence for Null 1.000 42993 29883
Daily individual’s experienced pushing 0.04 0.04 [-0.03, 0.11] 0.871 [-0.11, 0.11] 0.969 0.010 Very Strong Evidence for Null 1.000 51486 34192
Daily partner’s experienced pushing 0.11* 0.04 [ 0.03, 0.18] 0.994 [-0.11, 0.11] 0.527 0.177 Moderate Evidence for Null 1.000 38703 29065
Day 0.23** 0.08 [ 0.08, 0.38] 0.999 [-0.11, 0.11] 0.053 0.461 Weak Evidence for Null 1.000 71464 27970
Own action plan NA NA NA NA NA NA NA NA NA NA NA
Partner action plan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers -0.22*** 0.04 [-0.29, -0.16] 1.000 [-0.11, 0.11] 0.001 >100 Overwhelming Evidence 1.000 82288 27554
Daily facilitators 0.11*** 0.03 [ 0.05, 0.17] 1.000 [-0.11, 0.11] 0.524 2.455 Weak Evidence 1.000 79511 28190
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.13 0.29 [-0.45, 0.71] 0.671 [-0.11, 0.11] 0.273 0.015 Very Strong Evidence for Null 1.001 7580 14146
Mean partner’s experienced persuasion 0.42 0.29 [-0.15, 1.01] 0.928 [-0.11, 0.11] 0.106 0.042 Strong Evidence for Null 1.001 7465 13532
Mean individual’s experienced pressure -0.14 0.30 [-0.73, 0.45] 0.681 [-0.11, 0.11] 0.259 0.010 Very Strong Evidence for Null 1.001 9421 16495
Mean partner’s experienced pressure -0.31 0.29 [-0.90, 0.27] 0.855 [-0.11, 0.11] 0.170 0.016 Very Strong Evidence for Null 1.001 9111 16346
Mean individual’s experienced pushing 0.23 0.42 [-0.59, 1.05] 0.707 [-0.11, 0.11] 0.183 0.014 Very Strong Evidence for Null 1.000 9470 16049
Mean partner’s experienced pushing -0.06 0.42 [-0.88, 0.77] 0.553 [-0.11, 0.11] 0.205 0.012 Very Strong Evidence for Null 1.000 9488 16976
Mean barriers -0.56*** 0.08 [-0.72, -0.40] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 73715 31846
Mean facilitators 0.31*** 0.06 [ 0.19, 0.43] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 43106 32860
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.08 [0.48, 0.79] NA NA NA NA NA 1.000 11012 17933
sd(Daily individual’s experienced persuasion) 0.04 0.03 [0.00, 0.12] NA NA NA NA NA 1.000 12691 19187
sd(Daily partner’s experienced persuasion) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 15038 19782
sd(Daily individual’s experienced pressure) 0.07 0.06 [0.00, 0.27] NA NA NA NA NA 1.000 20707 22132
sd(Daily partner’s experienced pressure) 0.09 0.07 [0.00, 0.28] NA NA NA NA NA 1.000 18565 19559
sd(Daily individual’s experienced pushing) 0.07 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 13505 16761
sd(Daily partner’s experienced pushing) 0.10 0.05 [0.01, 0.21] NA NA NA NA NA 1.000 12544 11277
Additional Parameters
sigma 0.87 0.02 [0.84, 0.90] NA NA NA NA NA 1.000 61624 29562 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.73), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.73), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.76), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.84). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb +
    facilitators_self_cw + facilitators_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##       pressure_self_cw 3.96 [3.76, 4.19]         1.99      0.25
##    pressure_partner_cw 1.43 [1.37, 1.49]         1.19      0.70
##        pushing_self_cw 1.34 [1.30, 1.40]         1.16      0.74
##     pushing_partner_cw 1.17 [1.13, 1.21]         1.08      0.86
##     persuasion_self_cb 1.11 [1.07, 1.15]         1.05      0.90
##  persuasion_partner_cb 1.05 [1.03, 1.10]         1.03      0.95
##       pressure_self_cb 1.35 [1.31, 1.41]         1.16      0.74
##    pressure_partner_cb 1.18 [1.14, 1.23]         1.09      0.85
##        pushing_self_cb 4.64 [4.39, 4.90]         2.15      0.22
##       barriers_self_cb 4.92 [4.66, 5.20]         2.22      0.20
##   facilitators_self_cw 2.07 [1.98, 2.17]         1.44      0.48
##   facilitators_self_cb 2.25 [2.15, 2.36]         1.50      0.44
##                    day 1.33 [1.28, 1.38]         1.15      0.75
##  Tolerance 95% CI
##      [0.24, 0.27]
##      [0.67, 0.73]
##      [0.71, 0.77]
##      [0.82, 0.89]
##      [0.87, 0.93]
##      [0.91, 0.97]
##      [0.71, 0.77]
##      [0.82, 0.88]
##      [0.20, 0.23]
##      [0.19, 0.21]
##      [0.46, 0.51]
##      [0.42, 0.47]
##      [0.72, 0.78]
## 
## Moderate Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 7.80 [7.36, 8.26]         2.79      0.13
##  persuasion_partner_cw 9.24 [8.71, 9.79]         3.04      0.11
##     pushing_partner_cb 5.69 [5.38, 6.02]         2.38      0.18
##       barriers_self_cw 5.18 [4.90, 5.48]         2.28      0.19
##  Tolerance 95% CI
##      [0.12, 0.14]
##      [0.10, 0.11]
##      [0.17, 0.19]
##      [0.18, 0.20]
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'print': Predictive errors are not defined for ordinal or categorical models.
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 486, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.002057613
## sample estimates:
## outlier frequency (expected: 0.000432098765432099 ) 
##                                                   0
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 
## Sampling priors, please wait...
summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.99*** 1.46 [ 1.99, 8.38] 1.000 [0.84, 1.20] 0.000 52.547 Very Strong Evidence 1.000 28599 29417
Intercept[2] 8.62*** 3.28 [ 4.21, 18.79] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 28573 29205
Intercept[3] 23.84*** 9.75 [ 10.98, 55.70] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 27357 30723
Intercept[4] 100.41*** 48.37 [ 40.83, 271.18] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 29405 29297
Intercept[5] 1727.54*** 1494.86 [383.42, 12639.58] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 43980 30494
Within-Person Effects
Daily individual’s experienced persuasion 0.74* 0.09 [ 0.57, 0.93] 0.993 [0.84, 1.20] 0.149 1.611 Weak Evidence 1.000 29331 27721
Daily partner’s experienced persuasion 1.02 0.13 [ 0.77, 1.30] 0.558 [0.84, 1.20] 0.834 0.059 Strong Evidence for Null 1.000 37449 26256
Daily individual’s experienced pressure 1.78* 0.39 [ 1.09, 2.79] 0.986 [0.84, 1.20] 0.044 1.155 Weak Evidence 1.000 28334 24089
Daily partner’s experienced pressure 1.24 0.39 [ 0.53, 2.41] 0.742 [0.84, 1.20] 0.330 0.098 Strong Evidence for Null 1.000 21622 15563
Daily individual’s experienced pushing 1.31* 0.17 [ 1.02, 1.69] 0.981 [0.84, 1.20] 0.231 0.561 Weak Evidence for Null 1.000 39000 30661
Daily partner’s experienced pushing 0.90 0.13 [ 0.66, 1.22] 0.755 [0.84, 1.20] 0.667 0.072 Strong Evidence for Null 1.000 37897 28474
Day 1.40 0.67 [ 0.55, 3.56] 0.761 [0.84, 1.20] 0.231 0.066 Strong Evidence for Null 1.000 46914 30902
Own action plan NA NA NA NA NA NA NA NA NA NA NA
Partner action plan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.85 0.18 [ 0.55, 1.29] 0.774 [0.84, 1.20] 0.469 0.066 Strong Evidence for Null 1.000 52827 30927
Daily facilitators 0.82 0.17 [ 0.54, 1.22] 0.839 [0.84, 1.20] 0.426 0.091 Strong Evidence for Null 1.000 45213 30028
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.31 0.95 [ 0.31, 5.64] 0.647 [0.84, 1.20] 0.185 0.081 Strong Evidence for Null 1.000 22219 26063
Mean partner’s experienced persuasion 2.36 1.87 [ 0.52, 12.75] 0.867 [0.84, 1.20] 0.104 0.132 Moderate Evidence for Null 1.000 21615 23875
Mean individual’s experienced pressure 2.71 2.38 [ 0.48, 16.28] 0.875 [0.84, 1.20] 0.085 0.125 Moderate Evidence for Null 1.000 26090 29259
Mean partner’s experienced pressure 0.83 0.69 [ 0.15, 3.98] 0.591 [0.84, 1.20] 0.172 0.064 Strong Evidence for Null 1.000 22181 27225
Mean individual’s experienced pushing 0.76 0.77 [ 0.10, 6.21] 0.605 [0.84, 1.20] 0.135 0.072 Strong Evidence for Null 1.000 20225 27189
Mean partner’s experienced pushing 0.04* 0.05 [ 0.00, 0.49] 0.994 [0.84, 1.20] 0.005 2.048 Weak Evidence 1.000 22454 27000
Mean barriers 1.14 0.72 [ 0.33, 3.97] 0.584 [0.84, 1.20] 0.223 0.067 Strong Evidence for Null 1.000 28199 28866
Mean facilitators 0.89 0.33 [ 0.43, 1.87] 0.631 [0.84, 1.20] 0.355 0.082 Strong Evidence for Null 1.000 28117 27832
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.11 0.29 [0.63, 1.81] NA NA NA NA NA 1.000 12805 22280
sd(Daily individual’s experienced persuasion) 0.34 0.19 [0.02, 0.71] NA NA NA NA NA 1.001 5885 9410
sd(Daily partner’s experienced persuasion) 0.19 0.15 [0.01, 0.53] NA NA NA NA NA 1.000 15675 18239
sd(Daily individual’s experienced pressure) 0.46 0.32 [0.03, 1.26] NA NA NA NA NA 1.000 11216 14696
sd(Daily partner’s experienced pressure) 0.63 0.52 [0.03, 2.04] NA NA NA NA NA 1.001 9533 17764
sd(Daily individual’s experienced pushing) 0.22 0.17 [0.01, 0.63] NA NA NA NA NA 1.000 10476 18263
sd(Daily partner’s experienced pushing) 0.17 0.15 [0.01, 0.65] NA NA NA NA NA 1.000 17755 20914
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and
##   b_Intercept[2] (r = 0.79), b_Intercept[4] and b_Intercept[3] (r =
##   0.85), b_pressure_partner_cb and b_persuasion_partner_cb (r =
##   0.74). This might lead to inappropriate results. See 'Details' in
##   '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="06_SensitivityBarriersFacilitators_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.10 [1.07, 1.14]         1.05      0.91
##  persuasion_partner_cw 1.22 [1.18, 1.27]         1.10      0.82
##       pressure_self_cw 1.04 [1.01, 1.09]         1.02      0.96
##    pressure_partner_cw 1.05 [1.02, 1.10]         1.02      0.95
##        pushing_self_cw 1.10 [1.07, 1.15]         1.05      0.91
##     pushing_partner_cw 1.22 [1.18, 1.27]         1.10      0.82
##     persuasion_self_cb 3.13 [2.97, 3.31]         1.77      0.32
##  persuasion_partner_cb 3.46 [3.28, 3.66]         1.86      0.29
##       pressure_self_cb 2.16 [2.06, 2.27]         1.47      0.46
##    pressure_partner_cb 2.08 [1.98, 2.18]         1.44      0.48
##        pushing_self_cb 2.47 [2.35, 2.60]         1.57      0.40
##     pushing_partner_cb 2.33 [2.22, 2.45]         1.53      0.43
##       barriers_self_cw 1.31 [1.26, 1.36]         1.14      0.76
##       barriers_self_cb 1.25 [1.20, 1.30]         1.12      0.80
##   facilitators_self_cw 1.27 [1.22, 1.32]         1.13      0.79
##   facilitators_self_cb 1.08 [1.05, 1.13]         1.04      0.93
##                    day 1.07 [1.04, 1.11]         1.03      0.94
##  Tolerance 95% CI
##      [0.88, 0.94]
##      [0.79, 0.85]
##      [0.92, 0.99]
##      [0.91, 0.98]
##      [0.87, 0.93]
##      [0.79, 0.85]
##      [0.30, 0.34]
##      [0.27, 0.30]
##      [0.44, 0.49]
##      [0.46, 0.50]
##      [0.38, 0.42]
##      [0.41, 0.45]
##      [0.73, 0.79]
##      [0.77, 0.83]
##      [0.76, 0.82]
##      [0.89, 0.95]
##      [0.90, 0.96]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         34%
##          beta         16%
##        cauchy         16%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         75%
##  beta-binomial         16%
##       binomial          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 34 observations with a pareto_k > 0.7 in model 'model'.
## We recommend to set 'moment_match = TRUE' in order to perform moment
## matching for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with
##  approximate expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 486, p-value =
## 0.6217
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000520929 0.0114105115
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.002057613
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summary_is_reactance <- summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.55 0.22 [0.25, 1.21] 0.933 [0.83, 1.20] 0.119 0.221 Moderate Evidence for Null 1.000 30950 31371
Within-Person Effects
Daily individual’s experienced persuasion 0.70* 0.11 [0.49, 0.94] 0.989 [0.83, 1.20] 0.119 1.324 Weak Evidence 1.000 27108 26496
Daily partner’s experienced persuasion 1.13 0.20 [0.80, 1.70] 0.760 [0.83, 1.20] 0.585 0.102 Moderate Evidence for Null 1.000 31167 28174
Daily individual’s experienced pressure 2.32* 1.08 [1.03, 8.64] 0.978 [0.83, 1.20] 0.042 1.098 Weak Evidence 1.000 20159 20217
Daily partner’s experienced pressure 1.84 1.34 [0.44, 14.77] 0.822 [0.83, 1.20] 0.140 0.235 Moderate Evidence for Null 1.000 19691 18541
Daily individual’s experienced pushing 1.43* 0.23 [1.07, 2.08] 0.991 [0.83, 1.20] 0.120 1.212 Weak Evidence 1.000 29516 29930
Daily partner’s experienced pushing 0.89 0.22 [0.54, 1.53] 0.678 [0.83, 1.20] 0.493 0.106 Moderate Evidence for Null 1.000 32280 27957
Day 1.30 0.70 [0.45, 3.70] 0.682 [0.83, 1.20] 0.234 0.066 Strong Evidence for Null 1.000 49568 31621
Own action plan NA NA NA NA NA NA NA NA NA NA NA
Partner action plan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.76 0.19 [0.46, 1.25] 0.857 [0.83, 1.20] 0.327 0.101 Moderate Evidence for Null 1.000 52244 31167
Daily facilitators 0.78 0.19 [0.49, 1.25] 0.847 [0.83, 1.20] 0.364 0.107 Moderate Evidence for Null 1.000 47265 31320
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 3.45 3.60 [0.44, 30.97] 0.883 [0.83, 1.20] 0.068 0.226 Moderate Evidence for Null 1.000 22173 26252
Mean partner’s experienced persuasion 6.29 7.10 [0.73, 73.02] 0.953 [0.83, 1.20] 0.033 0.419 Weak Evidence for Null 1.000 20866 25481
Mean individual’s experienced pressure 56.19* 108.77 [1.65, 3647.05] 0.988 [0.83, 1.20] 0.006 1.592 Weak Evidence 1.001 19843 26744
Mean partner’s experienced pressure 0.52 1.12 [0.01, 32.15] 0.617 [0.83, 1.20] 0.062 0.171 Moderate Evidence for Null 1.000 15177 22911
Mean individual’s experienced pushing 0.48 0.81 [0.02, 13.35] 0.670 [0.83, 1.20] 0.078 0.126 Moderate Evidence for Null 1.000 17977 24375
Mean partner’s experienced pushing 0.03 0.05 [0.00, 1.02] 0.974 [0.83, 1.20] 0.011 0.809 Weak Evidence for Null 1.000 17773 26191
Mean barriers 0.71 0.60 [0.12, 3.59] 0.661 [0.83, 1.20] 0.156 0.096 Strong Evidence for Null 1.000 25812 28325
Mean facilitators 0.93 0.47 [0.34, 2.60] 0.558 [0.83, 1.20] 0.278 0.109 Moderate Evidence for Null 1.000 25527 27451
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 2.01 0.43 [1.30, 3.04] NA NA NA NA NA 1.000 14792 23932
sd(Daily individual’s experienced persuasion) 0.47 0.22 [0.06, 0.95] NA NA NA NA NA 1.001 8698 10543
sd(Daily partner’s experienced persuasion) 0.42 0.25 [0.04, 1.00] NA NA NA NA NA 1.000 14385 16368
sd(Daily individual’s experienced pressure) 1.16 0.68 [0.12, 2.90] NA NA NA NA NA 1.001 9847 10596
sd(Daily partner’s experienced pressure) 1.68 1.05 [0.14, 4.23] NA NA NA NA NA 1.000 14248 13579
sd(Daily individual’s experienced pushing) 0.29 0.23 [0.02, 0.83] NA NA NA NA NA 1.000 10535 17156
sd(Daily partner’s experienced pushing) 0.44 0.33 [0.03, 1.28] NA NA NA NA NA 1.000 13637 16986
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="06_SensitivityBarriersFacilitators_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

hypothesis(is_reactance, "exp(pressure_self_cw) > exp(pushing_self_cw)")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper
## 1 (exp(pressure_sel... > 0     1.47      2.38    -0.43     5.15
##   Evid.Ratio Post.Prob Star
## 1       5.45      0.84     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = FALSE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
Hurdle Component
Non-Zero Component
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 2.50*** [ 1.51, 4.14] 49.95*** [42.68, 58.63] 122.68*** [108.65, 138.51] 3.80*** [ 3.59, 4.02] 0.55 [0.25, 1.21]
Within-Person Effects
Daily individual’s experienced persuasion 1.46*** [ 1.23, 1.79] 1.02 [ 0.97, 1.09] 1.02 [ 0.98, 1.06] 0.01 [-0.04, 0.06] 0.70* [0.49, 0.94]
Daily partner’s experienced persuasion 1.32** [ 1.10, 1.64] 1.02 [ 0.98, 1.07] 1.03 [ 0.99, 1.07] 0.01 [-0.04, 0.06] 1.13 [0.80, 1.70]
Daily individual’s experienced pressure 0.80 [ 0.47, 1.49] 0.88 [ 0.77, 1.01] 0.95 [ 0.86, 1.05] -0.05 [-0.18, 0.08] 2.32* [1.03, 8.64]
Daily partner’s experienced pressure 1.52 [ 0.82, 5.09] 0.95 [ 0.86, 1.04] 0.96 [ 0.89, 1.04] -0.06 [-0.19, 0.06] 1.84 [0.44, 14.77]
Daily individual’s experienced pushing 1.17 [ 0.88, 1.60] 0.97 [ 0.91, 1.03] 1.03 [ 0.97, 1.09] 0.04 [-0.03, 0.11] 1.43* [1.07, 2.08]
Daily partner’s experienced pushing 1.60** [ 1.18, 2.47] 0.97 [ 0.91, 1.02] 1.01 [ 0.96, 1.06] 0.11* [ 0.03, 0.18] 0.89 [0.54, 1.53]
Day 0.97 [ 0.63, 1.50] 1.00 [ 0.88, 1.14] 0.97 [ 0.88, 1.07] 0.23** [ 0.08, 0.38] 1.30 [0.45, 3.70]
Own action plan NA NA NA NA NA NA NA NA NA NA
Partner action plan NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.71*** [ 0.58, 0.86] 1.22*** [ 1.13, 1.31] 0.99 [ 0.94, 1.03] -0.22*** [-0.29, -0.16] 0.76 [0.46, 1.25]
Daily facilitators 1.31** [ 1.07, 1.61] 1.23*** [ 1.17, 1.29] 1.06** [ 1.02, 1.10] 0.11*** [ 0.05, 0.17] 0.78 [0.49, 1.25]
Daily wear time NA NA NA NA 1.00*** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.83 [ 0.28, 2.44] 1.06 [ 0.77, 1.45] 1.01 [ 0.73, 1.39] 0.13 [-0.45, 0.71] 3.45 [0.44, 30.97]
Mean partner’s experienced persuasion 1.27 [ 0.43, 3.69] 0.92 [ 0.67, 1.26] 0.88 [ 0.63, 1.22] 0.42 [-0.15, 1.01] 6.29 [0.73, 73.02]
Mean individual’s experienced pressure 0.32 [ 0.07, 1.28] 0.65 [ 0.30, 1.42] 1.08 [ 0.75, 1.55] -0.14 [-0.73, 0.45] 56.19* [1.65, 3647.05]
Mean partner’s experienced pressure 0.25* [ 0.05, 1.00] 0.49 [ 0.23, 1.09] 1.03 [ 0.72, 1.46] -0.31 [-0.90, 0.27] 0.52 [0.01, 32.15]
Mean individual’s experienced pushing 1.09 [ 0.22, 5.34] 1.48 [ 0.83, 2.64] 1.10 [ 0.69, 1.77] 0.23 [-0.59, 1.05] 0.48 [0.02, 13.35]
Mean partner’s experienced pushing 1.48 [ 0.30, 7.25] 1.70 [ 0.95, 3.05] 1.29 [ 0.81, 2.07] -0.06 [-0.88, 0.77] 0.03 [0.00, 1.02]
Mean barriers 0.54** [ 0.34, 0.86] 0.94 [ 0.80, 1.09] 0.79*** [ 0.71, 0.88] -0.56*** [-0.72, -0.40] 0.71 [0.12, 3.59]
Mean facilitators 1.28 [ 0.91, 1.79] 1.10 [ 0.99, 1.21] 0.95 [ 0.88, 1.03] 0.31*** [ 0.19, 0.43] 0.93 [0.34, 2.60]
Mean wear time NA NA NA NA 1.00* [ 1.00, 1.00] NA NA NA NA
Random Effects
sd(Intercept) 1.23 [0.92, 1.67] 0.29 [0.22, 0.40] 0.33 [0.26, 0.44] 0.61 [0.48, 0.79] 2.01 [1.30, 3.04]
sd(Daily individual’s experienced persuasion) 0.16 [0.01, 0.43] 0.11 [0.07, 0.17] 0.05 [0.02, 0.09] 0.04 [0.00, 0.12] 0.47 [0.06, 0.95]
sd(Daily partner’s experienced persuasion) 0.22 [0.01, 0.52] 0.06 [0.01, 0.12] 0.05 [0.00, 0.10] 0.04 [0.00, 0.11] 0.42 [0.04, 1.00]
sd(Daily individual’s experienced pressure) 0.44 [0.02, 1.84] 0.09 [0.00, 0.33] 0.06 [0.00, 0.23] 0.07 [0.00, 0.27] 1.16 [0.12, 2.90]
sd(Daily partner’s experienced pressure) 0.82 [0.05, 2.38] 0.05 [0.00, 0.18] 0.04 [0.00, 0.13] 0.09 [0.00, 0.28] 1.68 [0.14, 4.23]
sd(Daily individual’s experienced pushing) 0.42 [0.04, 0.94] 0.06 [0.00, 0.15] 0.09 [0.02, 0.16] 0.07 [0.00, 0.18] 0.29 [0.02, 0.83]
sd(Daily partner’s experienced pushing) 0.42 [0.03, 1.15] 0.05 [0.00, 0.13] 0.05 [0.00, 0.12] 0.10 [0.01, 0.21] 0.44 [0.03, 1.28]
Additional Parameters
sigma NA NA 0.66 [0.63, 0.68] 0.54 [0.52, 0.56] 0.87 [0.84, 0.90] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()
  • rmarkdown (version 2.29; Allaire J et al., 2024)
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • rlang (version 1.1.4; Henry L, Wickham H, 2024)
  • tidybayes (version 3.0.7; Kay M, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.2; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • plotly (version 4.10.4; Sievert C, 2020)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()